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Introduction 


The major focus of this research work has been the continuing develop- 
ment of "Riccatl Iteration." a new numerical tool for the design end analysis 
of linear control systems. This report Is arranged Into four chapters con- 
sisting of two research papers, an extensive literature review, and a 
viewgraph summary of directions for future research work. The Masters degree 
thesis "Comparison of Algebraic Riccatl Equation solvers" by Mr. Rasim 
Baykan was also supported by tNs research grant. An apology must be made 
for the fact that the numerical examples In chapters II and III are taken 
from other aerospace applications and computer-generated random cases rather 
than aerodynamic flutter models as originally proposed In this research effort. 


POLE PLACEMENT AND ORDER REDUCTION IN TMO-TIME-SCALE 
CONTROL SYSTEMS THROUGH RICCATI ITERATION 

Leonard R. Anderson* 


Abstract 

A trinsformatlon of variables taken from singular perturbations may 
be applied to two*t1me-sca1e linear systems in state space form to reduce 
the system to block-diagonal form with slow and fast modes decoupled. The 
transformation Is easily computed by applying the new "Rlccatl Iteration." 

The Iteration yields a solution to the nonsymnetrlc algebraic Rlccatl equation 
obtained by partitioning the original system matrix A. The numerical pro- 
cedure Is Initiated with the trivial Iterate Lo * 0, and Is globally convergent 
to the desired unique time scale decoupling solution. 

After transfoiinatlon, *hr decoupled system may be used in controller 
design to achieve exact closed-loop pole placement In the slow subsystem 
without altering the poles of the fast subsystem. The decoupled form may 
also be used to reduce system order by sitting a small parameter to zero. 
Provided the fast subsystem Is stable* the order reduction can be expected 
to yield a good approximation to the original system. These methods are 
demonstrated using the 16th order linear model of a turbofan engine. 


This paper was presented at the 3rd International Conference on Mathematical 
Conference, University of Southern California, July 29-31, 1981 and has been 
accepted for publication In Mathematical Modeling. 


* Assist. Prof., Aerospace and Ocean Eng. Dept., VPI & SU 
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1 . INTRODUCTION 


Large scale linear control systems are encountered frequently In engineering 
problems such as the design of high performance a 1 rcraftU 2 , laroe space 
structures3, or power systems^, A recurring theme In such probl^s Is the 
need to reduce the originally high order of the system to facilitate simula- 
tion, elgenanalysis, or control system design studies. The question of order 
reduction has been widely researched and many order reduction schemes have 
been proposed5*D. The method presented here Is related to the familiar eigen- 
space analysis of linear multivariable systems, but Is particularly well suited 
to direct numerical Implementation for large scale systems, and systems with 
“stiff” dynamics. 

Consider the linear control system 


X «* Ax + Bu 

(1) 

w ■ Cx + Du 

where x, u and w are state, control and output variables of dimension n, m 
and 1, respectively. As proposed by Chow and Kokotovlc^, system (1) will be 
classified as two-time-scale If the eigenvalue spectrum of the A matrix, 
represented as x(A), can be separated by absolute values Into nonempty sets 
S and F with n^ and n 2 ^ n - n] elements, respectively, such that 

|s^| « |fjl for all s^ In S, and fj In F. (2) 

A naturally occurring system small paramete r Is the eigenvalue ratio 


y 


max|s^| 

m7F[7T"" ^ 

j ^ 


(3) 


This parameter provides a measure of the system's time scale separation and 
identifies (1) as a singular perturbations problem. This particular small 
parameter was also proposed by KelleyS. 


More refined partitioning of the spectrum A(A) yielding three or more tinve 
scales (i.e. • eigenvalue groups) may also be considered with two or more 
small parameters analogous to (3). The technique described below may be 
applied repeatedly for such cases. Only the two-time-scale case will be con- 
sidered here. 


2. THE LK TRANSFORMATION 


Such two-time-scale systems may be con'enlently transformed into decoupled 
subsystems by a two step transformatio t to new variables y * T 2 Ti x * Tx 
with 

n oi 


6 


r. 


K 


(5) 




I 


T - 


+ KL 
L 


K 


I 


( 6 ) 


where the I ar«> Identity matrices of dimension n<. and n2* If the original 
system (1) Is partitioned as 


^1 

• 

s 

'^1 

A, 2 

h 

+ 

8l' 

.*2. 


.^1 

Ag2^ 

X2 


Bz 

• m 


(7) 


where Is n] x n^, etc., then the new system defined by transformation 
(4) will be block-triangular provided the n2 x ni matrix L satisfies the 
nonsymnetrlc algebraic Riccatl equation (ARE) 


LAi 1 “ ^2^ ” ^12^ ^ ^21 * 

If, In addition the n«| x n2 matrix K satisfies the Lyapunov equation 

KA 22 - AiiK + A^2 * (9) 

where An * An - A12L, A22 * A22 + LA] 2* then ( 7 ) Is transformed to the 
block-diagonal form 


• • « 

^1 

• 

9 

" a ,, 0 



+ 

'V 

^2 


1 

0 

rT 


/2. 


h 

m m 


where = (I+KL) 8^+1(82*82“ L8^ + 82. 

Such transformations have appeared In singular perturbations work®*^^*^^ and 
were Introduced Into the controls literature by Kokotov1c^2. One attractive 
feature of transformation T Is that It Is always nonsingular and has the 
explicit inverse 



-K • 
I+LK 


Singe A Is similar to matrix A, the eigenvalues of A become the eigenvalues 
of A^ ^ and ^22* 1*e» 

x(A) “ x(A^^) U \{^ 22 ^ ’ 

13 14 

It is well known that ARE has many solutions If we express the A 

matrix in spectral form as A ■ MJQ where M Is a modal matrix, J Is diagonal 
(assuming A nondefective), and Q^M*', then solutions to ARE can be expressed 
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In terms of partitions of M. MJQ 1s partitioned conformably with (7) so that 
Ji and Qn have dimensions ni x ni, etc. It can be shown by matrix 
algebra that L Is a solution to ARE If 


1) Is nonsingular 

2 ) ^ ■ Q22 

Also, If L satisfies (11), then An - MnJiMn"'. In general there will be 
one solution to ARE corresponding to each partition of che eigenvalues of 
A Into A^^ and A 22 * 

The solution L of particular Interest In twootlme-scale systems Is the unique 
solution to ARE which yields 

x(A^^) • S, X(A22) ■ F (12) 

and from here on It Is assumed that L satisfies ARE gnd (12). The Lyapunov 
equation (9) will have a unique solutlonl® provided An and ^2 ^^®ve no common 
eigenvalues, which Is guaranteed by (12). The solution K of (9) can also be 
expressed In terms of the modal matrix M as 


K - M 11 Q 12 » “^12^22* 

One method of computing L Is to compute n-| eigenvectors 



slow elgenspace and apply (11). However, as recently shown, 
computed as the limit of the Riccatl Iterat1onl ° algorithm 


I spanning the 
3tl6, L can be 


h +1 * ^^2 * '■1^2^ ^ ^^1^1 ^1^ 


initialized with Lq « 0 which is globally convergent to the desired decoupling 
L. If we define the residual 


R^ « L^A^^ - A22L^ - L^A^2h ^ ^1 

liRinN 

then In the limit as i — ■*“ v.- 

Thus L is readily computed for strongly two-t1me>sca1e systems (i.e. those 
with very small u). 

3. ORDER REDUCTION 


The original system (1) can therefore be considered as the decoupled 
subsystems 


^1 

^2 


^1^1 ®1 “ 

(13a) 

^22^2 ^ ®2 ^ 

(13b) 
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which modal tha system slow and fast dynamic parts, raspactivaly. The 
appropriate Initial condition for (13a, b) Is y(0) ■ Tx(0). The decoupled 
form can yield a saving In simulation studies since only the fast subsystem 
{13b) need by Integrated with a small time step. If the orlolnal system 
(1) Is Integrated directly, then one must treat the entire nth order system 
as fast to obtain accurate numerical solutions. 

As noted above, order reduction Is often required In the study of large 
scale systems, and this may be conveniently done using decoupled subsystems 
(13a, b). Using the spectral norm. 


IIAzzll 




so we can write (13a, b) In the standard singularly perturbed form 


y, - A„y, + Bi .. 
u y^ - i^2 * ®2 “ 

A ^ A ^ 

where hz *^^ZZ 83 ■ w Bg. The zeroth order approximation to (14) 

obtained by setting the small parameter to zero Is then 

* 

yi • Ajiy^ + 8, u 

H ■ '^2 ' ®2 “ 


(14) 


(15«) 

(15b) 


with Initial condition 


y^(0) * (I+KL)x^(0) + K X 2 ( 0 ) 


A 

and the approximation x 


to the original state x Is given by 


X 




r-K 



(15c) 


(15d) 


If the fast subsystem (13b) has eigenvalues with large negative real part, 
I.e. 


Re (fj) « - max I s^ I for all fj In F (16) 

then the n^th order approximation x can be expected to yield a good approxima- 
tion for X outside of an Initial boundary layer where fast, stable dynamics 
may predominate. 

Note that if one Is primarily interested In output variables w rather than 
state variables x, the reduced order model (15) may be written in a compact 
matrix form similar to (1). Transform and partition the output matrix as 

[C, C^] * C - CT’^ 
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to yield 


>1 • *11 >1 * >1 “ 
m w 

w ■ C-j + 0 u 


(17) 


Mhere 


* w e « 

0 * 0 " ^ 2^22 ^ 2 ' 


This order reduction method has been shown to compare favorably with other 
standard methodsl". 


4. POLE PUCEMENT 

Dynamic systems are often modeled as linear control systems for the purpose 
of designing a feedback controller to achieve specific c1osed*1oop eigenvalue 
locations. Through application of the LK transformation, an n*d1mens1ona1 
eigenvalue placement problem can be reduced to separate eigenvalue placement 
problems of dimension n<| and U 2 > 

It Is well known that If a constant linear control system Is controllable . 
I.e. 


rank (B AB ... * n, 


(18) 


then there exists at least one real m x n dimensional feedback matrix H such 
that the closed*loop eigenvalues given by 

a(a+bh) 


can be placed arbitrarily (as long as complex eigenvalues appear In conjugate 
pairs). If the original system (1) Is controllable. It can be shown by linear 
algebra that the slow and fast subsystems (13a. b) are also controllable. 
Transfer of observability from the original system to the decoupled subsystems 
follows by a dual argument. 


The design of a feedback matrix H to achieve a specific closed*loop eigen- 
value location then follows directly. Suppose It Is desired to relocate the 
n. slow open-loop eigenvalues S to n^ new eigenvalue locations S'. If one 
cAn find an m x n^ feedback matrix which satisfies 

(All +B^H^)*S'» (19) 


then system (13a, b) with feedback u » H^y^ becomes 


^1 

• 


/2. 



ijH, 


0 

Aj, 



^l' 


.^2. 


( 20 ) 


Thus the fast eigenvalues F are unchanged and the corresponding feedback 
matrix for (1) Is given by 


[H, 0] T. 


( 21 ) 
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If In addition to relocating the slow mode eigenvalues S It Is also necessary 
to relocate the fast open-loop eigenvalues F to closed-loop valuf^ F' then 
this can be done as follows. After solving for feedback matrix In (19) 
find the unique (assuming $' and F have no common eigenvalues) nz x n] matrix 
P satisfying the Lyapunov equation 

P(Ail+8lH^) - A22P + BgH^ ■ 0. 

The block- triangular system (20) can then be transformed to a new block- 
diagonal form 



where y 2 * Py] y 2 » ^11 * ^11 82 ■ PB^ B 2 and u ■ u - H^y^. 

Now relgcate fast eigenvalues F to new values F* by finding an m x n 2 feedback 
matrix H 2 satisfying 

x(A22 + B 2 H 2 ) * F'- (23) 

The appropriate gain matrix for the original system (1) Is then given by 

u * Hx. H- [(ii^ + H 2 P) H 2 ]t. (24) 

Thus, the block-decoupling transformation T can be used to exactly relocate 
both slow and fast eigenvalues via state feedback. 

Many variations of this approach are possible to facilitate the feedback matrix 
design task. For example. by repeated application of the block-decoupling 
transformation either of the decoupled subsystems (13a. b) could be further 
transformed Into block-diagonal form and a complete eigenvalue relocation could 
be achieved In three consecutive design steps analogous to the above two design 
steps (19) and (23). In many practical applications, the fast modes are all 
stable and well damped so only the slow eigenvalues need be relocated. 

5. A TURBOFAN EXAMPLE 

The example considered here Is the 16th order model of a turbofan engine which 
was the theme problem for a recent conference on control of linear multivariable 
systems^. This model Is the linearization of a detailed nonlinear simulation 
(at the sea level maxlimim non-afterburner thrust point). The state variables 
consist of shaft speeds, temperatures and pressures; the five control Inputs 
are fuel flow, nozzle area, two vane positions and compressor bleed; and the 
five output variables are net thrust, total airflow, a temperature and two 
stall margins. For ni • 15, 5 and 3, the resulting small parameters are 0.304, 
0.371 and 0.383. Since these tes are not particularly small relative to 
one, we might call this system weakly two- time-scale. 

Selecting n^ « 5, the L and K matrices were computed as described In Anderson^^. 
The response 0 ^ the full 16th order system and reduced 5th order system con- 
taining the slow dynamics are compared In Figures 1 and 2. The response of 


n 



toUl thrust and fan ^pacd to stop Inputs In futi flow rata and Initt guldt 
vant position danonstrata that ^od agroament Is achlavad bttwatn full and 
raducad ordar rasponsa axcapt during an Initial boundary layar translant 
whara fast dynasties ara significant. Tha rasponsa to a gulda vana stap In- 
put previoas a savara tast of tha raducad ordar modal sinca this control 
warlabla Is locatad at tha front of tha angina and soma tlma Is raquirad 
for Us af facts to propagate to tha nat thrust. 

Pcia placamant for this axampla was carried out with tha goal of Increasing 
the spaed of tha thrust rasponsa. Applying tha methods of tha previous 
section, tha three slowest modes -0.65, -1.90 and -2.62 were relocated to 
approximately -6.6. Since there ara five control Inputs In this axampla and 
B] Is nonsingular. It was posslbli to shift tha above three poles and leave 
tha eigenvectors of the slow subsystem unchanged. Let tha slow subsystem 
have spectral form 

A,, . H,,J, M,,-' 


Where • dlag (-0.65, -1.90, -2.62, -6.72 i ^ 1.31). 
Then feedback matrix was chosen as 


where c > dlag (-6, -5, -4, 0, 0) and the model response for the resulting 
H (c.f. equation (21)) is Illustrated In Figure 3. As expected, the response 
Is much qi^lcker with this feedback. 


6. CONCLUSIONS 

A method of time scale separation based on the block-decoupling LK trans- 
formation is proposed. Methods of applying this transformation to order 
reduction and pole placement design tasks are described and demonstrated 
for a turbofan engine example. As previously shown, the L and K matrices 
are easily computed for large scale systems, particularly those systems 
with large time scale separation and correspondingly small 
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Response of thrust and fan speed to a 500 Ib/hr step change In fuel flow rate 




Figure 
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Response of thrust and fan speed to a 10 deg. step change In Inlet guide vane position 




Figure 



15 


Open- and closed-loop response of tiirust and fan speed to step changes in fuel flow rate 
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NuRwrlcal Solution of the Symmetric Riccatl Equation through Riccatl Iteration 


Leonard R. Anderson* 

Dennis W. Breiner** 

A* Rasim Bay kan*** 

Abstract 

This research paper presents a new numerical method for solving the 
symmetric algebraic Riccatl equation from optimal control. This algorithm 
employs the "Riccatl Iteration" which has been successfully used to solve 
time-scale decoupling problems In structural vibrations. The algorithm Is 
related to the subspace Iteration method, and the rate of convergence to the 
solution Is governed by the relative separation between the stable and unstable 
eigenvalues In the Hamiltonian system of equations. Provided there Is adequate 
eigenvalue separation and Ignoring roundoff error, the algorithm Is globally 
convergent to the desired Riccatl solution. The method Is demonstrated for 
a set of 8th order random examples. Preliminary accuracy and timing comparisons 
with other standard methods of solving the symmetric Riccatl equations are 
presented. 


Submitted for presentation at the 1982 American Control Conference, Arlington, 
Virginia and for publication In the IEEE Transactions on Automatic Control . 
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♦♦♦Graduate Research Assistant, Aerospace and Ocean Eng. Dept., VPI & SU. 
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I. Introduction 

It Is mU knoMn that In order to solve the Infinite-time regulator 
problem [1,2] with linear system 

t ■ Ax + Bu (1) 

and quadratic performance Index (the LQR problem) 

J • 1/2 J^(x^Qx + u^Ru)dt (2) 

one must compute the unique positive definite solution S to the symmetric 
algebraic Riccatl equation 

0 - SA + A^S - SBR’U^S + 0 . (3) 

Here x Is an n-dimenslonal state vector, u Is an m-dimenslonal control vector, 

A and B are constant matrices Q and R are positive seml-definite and positive 
definite, respectively. Also It Is assumed that the matrix pair (A,B) Is 
controllable, and the pair (A,Q^^^) Is observable and that B has full rank. 

Then assuming that any constraints on the size of u(t) are not violated, the 
linear feedback control law which minimizes the performance Index J Is given 

-1 T 

u « Kx where K * -R B S . (4) 

Algebraic Riccatl equations similar to (3) also play a fundamental roll In 
discrete-time optimal control problems. In optional estimation problems, and 
Is many other areas of applied mathematics [3-7]. See Jones [8] for basic 
necessary and sufficient conditions for solutions of equations of this type. 

Among the most widely used solution methods for solving algebraic 
Riccatl equations are the Schur vector method described by Laub [9], methods 
employing eigenvectors of the Hamiltonian system [1,2,10,11] and Newton -Raphson 
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mthods C1|2,12]« The Riccatl Iteration method diffiers from these methods 
In that neither eigenvalue computation nor solution of equations of the 
Lyapunov- type are required. The algorithm Is suprisingly simple, requiring 
only repetitive solution of linear equations of order n. However . caution 
must be exercised when applying the method to LQR problems with nearly zero 
eigenvalues, or eigenvalues with large Imaginary part* 

II. The eigenvector solution to the algebraic Riccatl equation 

Given the linear system (1) and performance Index (2), we can form the 
control Hamiltonian and apply the necessary conditions for optimal control to 
obtain the Hamiltonian system of equations 


(5) 


R 


‘a -BR‘^8^ 


’ X' 


» * 
X 


m 




i H 


• 

y 


-Q -A^ 

1 « 


y 

m « 


7. 


where y 1s an n-dl mens Iona 1 vector of costate or adjoint variables. The 
2n X 2n matrix H has synnetric eigenvalues. That Is, If X Is an eigenvalue 
of H then -X Is also an eigenvalue of H. 

If matrix H Is represented In spectral form 


H - M A 


Mil ^12 


'■ -A, 0 ' 


^11 ^12 

M 21 M22^ 


0 A^ 

«■ 


Q 2 I Q22 


where the matrices are partitioned conformally with (5), and contains all 
the eigenvalues with positive real part ordered as 

A^ ■ dlag (X^, X 2 * .*., X^^) with |X^|^jX^^^|, 1 * 1,2, ... n-1, (6) 

then the desired solution, S, to the Riccatl equation Is given by 


(7) 


or Q 22 S « . 
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To damonstrate this fact, apply tha similarity transformation 



» m 

I 0 


0 

a 

II 





f 

a-a 

</l 

I 

a 


9 

C/» 


( 8 ) 


to tha Hamiltonian system to obtain 


H = THT 


-1 


A - BR'^B^S 


•sa-aVsbr‘U^s-q 


-BR*^B^ 


-A^+SBR"^B^ 


(9) 


or In spectral form 




”i2 


N. 


'21 


”22 


r * A , 0 


0 A, 


^12 


-^21 ^ 22 - 


1 - 022^21 ^ 




0 


”ll '^1^1 2*^1 2*^2^22’ 

Q22 Q22 


(10) 


Since the. closed- loop system for feedback law (4) Is 

X « Ax + 8u « (A-BR’^B^S)x (11) 

= 

the eigenvalues -A| become the eigenvalues X(A) of the closed-loop system. 

Given efficient elgenanalysis software such as EISPACK [13]« one can directly 
compute the eigenvalues and eigenvectors of the Hamiltonian system (5) and solve 
the appropriate system of linear equations (7) for the Riccatl solution S. 
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If •A'l Includes complex conjugate eigenvalues Xj, ^ with corresponding complex 
conjugate eigenvectors Uj ■ a -t* 16, Uj « a - 16, then complex arithmetic 
can be avoided in the set of linear equations (7) by replacing Uj with a and 
Uj with 6. 


As an alternative to computing eigenvectors 




, one can compute 


the n Schur vectors [9] spanning the stable elgenspace and solve (7) with 
M 2 ^ replaced by corresponding Schur vectors. Schur vectors are defined as 
follows. Given any nxn matrix C with eigenvalues x.|, ..., there exists a 

U 

unitary matrix V such that V CV Is upper triangular with diagonal elements 

X| X^. The columns of V are the Schur vectors of C, possibly complex, 

corresponding to this particular ordering of eigenvalues. 

If some of the eigenvalues X^ are complex, then one can avoid complex 
Schur vectors by employing a (real) orthogonal matrix V such that V AV Is 
real and nearly upper triangular. The only nonzero subdiagonal elements In 

Vr 

V AV will be due to real 2x2 diagonal blocks corresponding to complex eigenvalues 


'V< 

Xf. This Is termed the real Schur form of C, and the columns of V are real 
Schur vectors. 

To return to the problem of computing the solution (7) to the symmetric 
Riccatl equation (3), one can compute the2n real Schur vectors 





V » 


( 12 ) 



rv< 


of the 2nx2n matrix H so that 



provides a basis for the elgenspace 


corresponding to stable eigenvalues 

The computation Is best performed using the subroutine HQR3 [14] with 
small modifications. In Its standard form, HQR3 orders the eigenvalues of H 
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by modulus with |X^U|X 2 l< ... on the diagonel of The ordering statements 
In HQR3 may be modified so that either 

1) all stable (negative real part) eigenvalues are placed In the upper left 
diagonal block of and unstable eigenvalues are placed In the lower 
right block* with no regard to ordering within each block* or 

vr 

2) the eigenvalues of H are oroered along the diagonal of V HV by real 
part from most negative to most positive. 

Only modification (2) Is Implemented In the numerical examples which 
follow. 

III. Riccatl Iteration solution to the algebraic Riccatl equation 

Before proceeding to the solution of the symmetric Riccatl equation (3)* 
we will first describe In detail the solution of the closely related time- 
scale decoupling problon using Riccatl Iteration [15*16], and then extend 
these results to (3). 

be stated as follows. Given an nxn 


-l^ngl* 

with u » ' ■ t (13) 

I'll 

If parameter y Is small relative to one, we say that the matrix Is two-time- 
scale [17*18]. Let the diagonal Izatlon of C be given by C ■ EOF 
where F«E’^ and J 1$ diagonal with diagonal entries 
^1 » s^* »..* , f^» f2» ...f fog • 


The time-scale decoupling problem can 
diagonal Izable matrix C with eigenvalues 
s^ * s^* ... * s„ * f^* f^, ...*f^ 


1 


not necessarily distinct and satisfying 


|Sli*lS2l < ■■■ S|lh,l'‘l^lsl*'2ls 
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Partition the matrices C, Et P and J as 



■Cn Hz 


■^11 

^12 

c ■ 


. E • 




Cjl C 22 

m m 


^1 

1 

CM 

CM 

Ul 



■^11 

’'12' 


■Ji 

0 ■ 

F ■ 



. J ■ 




/21 

^22. 


0 

'’ 2 . 


where n.|Xn^, etc. For k * 0. 1, 2, ... partition the matrix 

i( 

C similarly as 


C 


k 





The Riccatl iteration and Its convergence propertlei^ for the time-scale 
decoupling problem [19] are introduced through 

Theorem 1 . Given w < 1 and Vj^ nonsingular for k » 1, 2, 3, .... the Ite'^tlon 


- 0 

'• k+1 ' ^‘'22 * '■ k '' 12 ' * *' 21 * 


(14) 


k « Ot 1, 2, ... 


Is well defined and 


•1 


^22 ^21 ^ 


with convergence In the sense of matrix norm. For any IxJ matrix Bt |jB|| 
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will be given by 


llBil ■ sup llBxllj 

iixili 

where ||*(l<|* ^^ote convenient norms on and R*^. 

The proof of Theorem 1 Is presented In the Appendix. 

Furthermore, the rate of convergence of L|^ to ^ 22’^ ^21 controlled 

by the small parameter y as shown In the Appendix. Convergence would be 
expected to be unreesona:^ 1 y slow for y close to one. 

To apply the Riccatl Iteration to the symmetric Riccatl equation (3) 
we need to shift the Hamiltonian matrix (5) as 

H ■ H + X 5 I (15) 


so that the stable, left half-plane eigenvalues -Xp -X 2 , -X^ of H 

become the slow or small eigenvalues of H analogous to Sp S 2 » s,^^. 

We can now state the formal conditions for convergence cf the Riccatl 
Iteration for the synnatrlc Riccatl equation (3) as 


Theorem 2. Given the Hamiltonian system (5) and given shift distance 
Xj satisfying^ 


1 


Note, this redefines y for the 


sequel, cf. (13). 
24 


u ■ 


MX 

1«1 n 

min 



thfi Riccatl Iteration 


(16) 


So • 0 


^ XjD-Q) 


(17) 


converges to solution (7) with asymptotic rate 
! |S^ + Q22 ^21 1 1 - 


where the matrix norm Is as In Theorem 1, and constant c Is Independent of k. 
The proof of Theorem 2 follows directly from Uos proof of Theorem 1 In the 
Hamiltonian setting. 

The key practical difficulty In Implementing Iteration (17) Is the 
estimation of the shift distance X . Two choices for the shift distance are 
proposed, namely 


X 1»X ix^l 

1«l,...,n 

or 


(18) 


Xj i max |X^ 1 + JD^n 

1*l,.t.,n 1^1, ...fO 


I 



(19) 
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Tht shifting strategy (19) will In general yield higher rates of convergence 
since the ftable eigenvalues of H become nearly centered about the origin 
In ?f. Shift strategy (18) Is recommended for symmetric Riccatl equations with 
eigenvalues having large Imaginary part, since this minimizes the chances of 
no convergence, or convergence to a solution other than the desired solution (7). 
Both shifting strategies can be Implemented numerically In approximate form. 

The computational algorithm for the symnetrlc Riccatl equation Is therefore 
Algorithm 1. 

1) Form the Hamiltonian matrix 
r A -BR'^bT 


H - 



2 ) Beginning with the 2n vector 
Xq - (2n)*^/^ Cl 1... 1]^ 

perform a few power steps (five or more) 


'\( T T 

^ 1+1 • *1 H 


^1+1 * *1+1^1 ^^1-fl J ^ ^ 


where IMI Is the Euclldsanor Frobenlus norm. 


3} (tttain an Initial estimate of the shift distance 


a, - ||x,ll 


H+1 


‘1 

ai ^ I 


,1*1 ,2 , . . . , 5 


^s,0 • °6* 
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4) Estimate a limit for 

^limit * hiO 

and partition the vector Xg as 

*6 * 

with of dimension n. 

5) Initialize Sq » 0, j » 0, and form the inverse power vector 

Vq « [1 l...lf 

of dimension n. 

6) Form the matrices 


«22,J ■ 


• Sj(A * - Q. 


7) Perform a power step and estimate 


ax 


“jli * hz,s 


W,j * ^^limit* ll“j+lll ’ 

Vi“ Viviii 

8) Solve the linear system for and Vj^^ 
”22,J t^j+1 ^J+1^ “ ’j’ 


( 20 ) 


9) If j > It find the relative change in Sj. 
«,,* IIV, -Sjll/IISjll 
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and tarmlnata the Iteration If 6j Is less than or equal to some specified 
convergence tolerance e. or has reached a minimum value. 

10) Normalize the Inverse power vector and estimate 
■ "»* (0 . min -x->— - 

Vi “ 


11) update the shift distance by either 


or 

, _ ®^s.j W.j \n1n.j 

"s.j+1 ^ 


( 21 ) 


( 22 ) 


12) Increase j by one and go to step 6. 

Obviously. If one computes the updated shift distance by the more 
conservative (21), then one need not Include the Inverse power vector Vj In 
the Iteration. The power and Inverse power estimates of « |X^ | and 
^tmx * based on the fact that the eigenvalue spectrum of H 2 , j converges 

to the eigenvalues Increased hy j. 

This Iterative technique for solving the symmetric Riccatl equation has 
similarities to one proposed by F«rrarand 01 Pietro [20]. The methods differ 
In that their method Is Initialized using eigenvectors, and employs Iterative 
Improvement requiring solution of a Lyapunov-type equation at each step. 

IV. Numerical examples 

The eigenvector and Schur vector solution methods and the Riccatl Iteration 
method were evaluated In a limited comparative test. The test consisted of 
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sets of 11near-quadrtt1c*reguUtor problems generated as follows. The elements 
of the A and B matrices of system (1) were chosen as random Integers from a 
uniform distribution over the Internal [-lOt '*'10] » with 2SS of the elements 
set arbitrarily to zero. The order n of the system was varied between 5 and 20* 
and the control dimension m was held at one. The positive semidefinite weighting 
matrix Q was constructed from a Choleski factorization 

Q ■ PP^ 

where P Is upper- triangular. The upper-triangular elements of P are also 
Integers chosen uniformly from the Internal [-10, +10], with 25X of these 
elements set to zero. In all cases the R matrix, a scalar, had a value of 2. 

All methods were either obtained or coded In the FORTRAN language In 
double precision, and executed In batch mode on the VPI&SU IBM 370/158 computer. 

Table 1 presents the results for a set of ten random cases of order eight 
generated as described above and solved by Algorithm 1. The table Includes 
the theoretical rate of convergence u for each of the shift strategies (18) 
and (19), the value of the Riccatl residual ||R^i| where 

R^ * S^A + a"^S^ - S^BR“^b”^S^ + Q (23) 

for the converged solution S^, the maximum relative error between the 
eigenvalues of the closed-loop system (11) and the stable eigenvalues of the 
Hamiltonian system (5), the number of Iterations In Algorithm 1 to reach a 
minimum value of 6^, and the corresponding central processor execution time 
In seconds. As noted earlier, the shift strategy (18), (21) Is significantly 
slower In conve>^gence than (19), (22). However, the faster method (19), (22) 
falls to converge for case 8 as expected since u > 1. 

To further describe the convergence characteristics of the Riccatl 
Iteration, the values of log^QKR^Il and log^Q||d^|| versus Iteration 


29 




number are plotted In Figures 1 and 2. These figures correspond to cases 2 
and 9, and represent the most rapid and the slowest convergence of the ten 
cases, respectively. As shown, the theoretical limit y, cf. (15), accurately 
predicts the true rate of convergence of the algorithm for each of the two 
shifting strategies. The specific A, B, Q and R matrices for cases 2 and 9 
are listed In Tables 2 and 3. 

For comparison purposes, the symmetric Riccatl equation for these ten 
cases Is also solved by the eigenvector and Schur vector methods. This data 
Is listed In Table 4. As In Table 1, we specify the value of the norm of the 
Riccatl residual |{R(| corresponding to the computed solution S, the maximum 
relative error between the eigenvalues of the closed-loop system and the 
Hamiltonian system, and the execution time In seconds. The execution times 
Include both the time to compute eigenvectors or Schur vectors and the time 
to solve the linear system (7). The numerical method used tc solve both the 
linear systems (7) and (20) In this study was LU decomposition with partial 
pivoting without Iterative Improvement of accuracy. Specifically, the Fortran 
subroutines DECOMP and SOLVE were used, but not the subroutine IMPRUV from the 
widely used Forsythe and Moler text [21]. 

For the eigenvector solution method, all eigenvalues and eigenvectors 
of the Hamiltonian matrix were computed using EISPACK subroutines ELMHES, 

ELTRAN and HQR2. Then eigenvectors corresponding to stable eigenvalues were 

Pill 

selected to form tte real matrix. For the Schur vector solution method, the 
‘y 1 ^21-1 

real Schur vectors ^ spanning the stable elgenspace were computed 

using EISPACK subroutines ORTHES and ORTRAN, and Stewart's subroutines HQR3 
and EXCHNG. The HQR3 subroutine was modified [22] so that the eigenvalues 
m are ordered along the diagonal from most negative to most positive real part. 

As shown In comparison of Tables 1 and 4, the Riccatl Iteration Is able to 
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produce final Riccatl residuals two orders of magnitude smaller on the average 
than the eigenvector and Schur vector methods. However, the Riccatl Iteration 
execution times are slower than the eigenvector and Schur vector methods by 
a factor of ten, on the average. 

The choice of method clearly depends on both the parameter \i for the class 
of problems under consideration, and the relative Importance of minimizing 
computer execution time versus software simplicity. In general, as the order 
n of the LQR problem Increases, the parameter u increases and approaches 
or exceeds one. Therefore, the elgenvector/Schur vector algorithms would 
still be the method of choice for problems of large order. 

V. Conclusions 

A new numerical algorithm for solving the symmetric Riccatl equation from 
the linear*quadratic*regulator problem has been presented and compared with 
standard methods. A formal proof of convergence for the Riccati iteration Is 
presented, and numerical examples confirm the theoretical rate of convergence. 

The strengths of this new algorithm are its simplicity, accuracy and theoretically 
transparent basis. Riccati iteration may be particularly useful for low order 
adaptive control algorithms or control system design studies where one must 
update the Riccati solution as system matrices or performance index weighting 
matrices slowly change. The primary weakness of this method is that it is 
slower than standard methods as shown here, and the rate of convergence is 
dependent on sufficient stable/unstable eigenvalue separation in the Hamiltonian 
system. 
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Future research tasks Include: 

- finding an appropriate conformal transformation for the Hamiltonian 
system to increase stable/ unstable eigenvalue separation and/or 
decrease the imaginary component of eigenvalues. 

- Applying a doubling algorithm [23] to this method to Increase rate 
of convergence. 

Appendix - Proof of Theorem 1 

To prove theorem 1 the following lemma will be useful. 

Leitina: For k = 0, 1 , 2, ... 


Proof : Ue proceed by Induction on k: 

For k a 0 

Lo • Vj’Uo . I.O = 0 

as required. Note that partitioning the product 
C*^C » 
yields 


‘^k^n ''k^Zl * \+l 

\^i2 V22 * Vr 

Assuming the Induction hypothesis, 

,-l 


(C22 * ^22 '^k'^k^l 2 * Vl 

which Is nonsingular by hypothesis. 


Therefore L 


k+1 


(C 22 + L,^C,2^ ^^^k^ll '*■ ^21^ 

<viliVk){v’k\cn * C2 t) 

''k+1 '^k'^k^^k+l “ ''k+lVl • 


This completes the proof of the lemma. 
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Proof of Theorem 1 

Since C * EOF. where F»E’^, we have ■ (EJF)*^ * 
Partitioning the product EJ*^F yields 


\ * ^ 21 ^ 22 '^ 2*^^21 
''k “ ^ 21 '^! *^^12 ^ 22 ‘^ 2*^^22 


where and J 2 are diagonal matrices with diagonal entries Sp $ 2 * **** 
and f^, f 2 » f„^ respectively. Therefore, by the lemma 

‘•k * ^^ 21 '^ 1*^^12 ^ 22 '^ 2 *'*^ 22 ^ ^^^ 21 '^ 1*^^11 


for k » 0, 1, 2, ... 
For brevity let 
D . 


^^ 21 '^ 1^2 ^ 22 '^ 2 *^’' 22 ^- 


An easy computation shows that 


therefore 


Lk ' (I - V^21 '^iS2^^^22'^2'^‘'22^^4i'^/''i1 * 

* (I * D|j ^21'^1 ^12^^^22'^2 ^22^21'^! ^11 ^ ^99^91)* 


22' 21 


Hence, 


‘•k ■ '' 22^21 * '"22'^2*^^2^21'^/''ll 


-(Dk^E2ij/Fl2)(p22‘^2 Sz^l'V Hl^ 
-( 0j ; lE2ij / Fi2 )^ f22 ^ 2l )* 


Note that (IJ^H < |s„^| and |jJ 2 ’^|| < j^j 




so 
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Furthermore, 






^lJlS2 


Now 11(1 +T)'^ II < (1-||T||)‘’ If ||T|| <1. 

SO 1 1 D|c^ ^1^/^1211 ^ 1^. ayk 

for sufficiently large k, 

where a - ||F 22 ’’|I llEjj-’ll IIE^,!! IIFjjH. 

Combining these estimates yields 

ll‘■k-’^22■’''^1N - l|F22''f2lll 

Where b » a| |F,^ | j/| |F^2l I • 

Since u < 1* 0 w k— ►». This estimate completes the proof. 
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Table 1 Summary of Riccatl Iteration Performance for ten 8th order random examples 




































































































Table 2 System and performance Index matrices for case 12 


A * 


Q 


‘ -8 

-8 

6 

-4 

-7 

4 

0 

0 ' 

-2 

-10 

-6 

-4 

0 

-3 

1 

1 

-5 

9 

2 

S 

-6 

8 

0 

-7 

-7 

-4 

-4 

10 

5 

2 

0 

-3 

3 

5 

-1 

2 

1 

-10 

0 

6 

-8 

2 

-7 

-9 

-8 

4 

-8 

0 

-7 

-1 

-6 

9 

-10 

-8 

-6 

7 

. -4 

7 

0 

4 

0 

6 

-1 

2 . 

C 0 

-7 

-7 

5 

6 

-7 

-2 

0.3 


■ 48 

18 

7 

-3 

-5 

-10 

6 

-9 ■ 

18 

27 

-10 

11 

-1 

3 

0 

0 

7 

-10 

30 

-6 

1 

-5 

0 

0 

-3 

11 

-6 

33 

0 

2 

0 

0 

-5 

-1 

1 

0 

6 

-3 

-4 

6 

-10 

3 

-5 

2 

-3 

10 

0 

0 

6 

0 

0 

0 

-4 

0 

4 

-6 

. -9 

0 

0 

0 

6 

0 

-6 

9 . 


R « 2 
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Table 3 System and Performance Index matrices for case #9 


A • 


Q - 


■ -8 

*8 

-7 

0 

•4 

10 

-9 0 ■ 

4 

-9 

-6 

-4 

7 

0 

-4 -7 

1 

7 

-10 

-1 

-3 

0 

10 -2 

0 

6 

10 

0 

0 

9 

2 -10 

0 

0 

0 

0 

0 

0 

3 0 

7 

-2 

-3 

2 

-4 

-6 

2 0 

8 

-8 

8 

0 

0 

0 

0 -8 

. 5 

0 

0 

0 

-2 

-1 

-5 -6 . 

CO 

5 

8 

6 

0 

0 

0 -6 3 


■ 44 

7 

-6 

15 

17 

-7 

-12 

6 ‘ 

7 

20 

8 

7 

7 

-11 

-4 

-2 

-6 

8 

21 

6 

-4 

-4 

4 

-4 

15 

7 

6 

29 

11 

-7 

-6 

4 

17 

7 

-4 

11 

14 

-7 

-10 

4 

-7 

-11 

-4 

-7 

-7 

17 

2 

4 

-12 

-4 

4 

-6 

-10 

2 

8 

-4 

6 

-2 

-4 

4 

4 

4 

-4 

4 . 


R « 2 
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Table 4 Swury of eigenvector and Schur vector solution perfonunce for ten 8th order randoa exa^>1es 




Fig. I The Convergence of Case 12 



log |R 
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A Literature Review of Robust Controller Design Methods 


Alok Oas* 


Abstract 

This paper presents a literature survey on the methods available for 
designing robust controllers. A mmber of methods for reducing the trajectory/ 
performance Index sensitivity In linear regulators are described. It Is shown 
with the help of an example that decrease In system sensitivity to variation 
In parameters Is obtained at the cost higher value for the performance 
Index. A method for reducing the eigenvalue sensitivity Is also discussed. 


* Graduate Research Assistant, Aerospace & Ocean Eng. Dept., VPI & SU 
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Introduction 


The design of a controller for a system requires a mathematical model of 
the system. In most cases some of the parameters of the system will change 
with time. The changes in the parameters could occur due to a variety of 
reasons for example, aging of components, environmental changes, etc. It has 
been found that for some systems even small changes In certain key parameters 
could appreciably degrade the performance of the control system. Hence, It 
is very important to estimate the effect of changes In the system parameters 
on its overall performance and to use controllers which minimize this degradation 
in performance. 

The problem then is, given a specified structure of the controller, to 
find the particular controller which yields a system with minimum sensitivity 
to variation In system parameters. Such controllers are usually called robust 
or Insensitive controllers. The robust controllers with constant feedback 
gain matrix will be considered. Adaptive controllers provide an attractive 
alternative to the problem of varying system parameters. We will not discuss 
adaptive controllers here. 

A large amount of literature Is available on robust controller design 

using both classical and modem control theories. A brief literature survey 

is presented here and some of the more Important methods are described. In 

the analysis, the system parameters are assumed to take an unknown but constant 

value around their known nominal value. 

Appendix A contains a bibliography on robust controllers containing 173 

references. The bibliography covers the following journals and conference 

proceedings from 19‘/0 to 1980. 

Allerton Conference on Circuit and System Theory 
IEEE Transactions on Automatic Control 
International Journal of Control 
Joint Automatic Control Conference 
IEEE Decision and Control Conference 

IFAC Symposium on System Stability and Adaptivity (1968, 1973) 

Automatica 
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Some other Important references are also Included. 


Parameter Sensitivity Reduction In Linear Regulators : 

A considerable amount of work has been done on the parameter sensitivity 
reduction In linear regulators. Most of the literature on this topic falls 
In one of the following categories: 

a) Trajectory sensitivity reduction 

b) Performance Index sensitivity reduction 

In both of these approaches to sensitivity reduction* the main objective 
Is to achieve a trade off between optimality In the nominal performance and 
sensitivity to small parameter variations. In a recent paper Yedavalll and 
Skelton [1] treat the problem of trajectory sensitivity and performance Index 
sensitivity In a unified way. We will now discuss the trajectory sensitivity 
reduction and combined trajectory and performance sensitivity reduction. 

Trajectory sensitivity reduction. In this approach a quadratic trajectory 
sensitivity term is Included In the Integrand of the performance index. One 
of the Initial papers on this Is by kahne [2]. This paper served as a starting 
point for a number of research efforts. Kahne Implements the control In an 
open loop manner. 

Let the linear time Invariant system be given by: 
x(t) » Ax(t) + Bu(t), x(o) = Xq (1) 

where x(t) is the state vector of dimension n 
u(t) Is the control vector of dimension m 

A and B are the state and control matrices of appropriate dimensions, 
a Is a time invariant parameter of the system x« A may depend on a. B 
and u are taken to be Independent of a. 

In the linear regulator problem we deten<i1ne the optimal control vector 
u (t) which minimizes the performance Index (subject to system eq. (1)}. 
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J 


1 

I 


x^(T) Fx(T) + 



Cx^(t) Q(t) x(t) + J(t) R(t) u(t)] dt 


( 2 ) 

( 2 ) 


where T Is the final time. 

F and Q(t) are positive semidefinite nxn matrices 
and R(t) is positive definite mxm matrix. 

To incorporate the trajectory sensitivity into the design, consider the 
trajectory sensitivity vector 


a * 


aa 


a s a nominal 


Differentiating (1) with respect to a gives 


a » Aa + A x, a(o) * 0 
a ' ' 


(3) 


where A * 

a 9a 


A trajectory sensitivity term is incorporated in the performance index 
and the problem becomes: 

find u(t) that minimizes 


J = ^ x^(T)Fx(T) + ~ |^{x'^Qx + u^Ru + dt (4) 

o 

subject to 


X = Ax + Bu, x(o) = Xq 
0 » Aa + A^x* o(o) = 0 

Kahne solved this problem in the standard way. 

This method has two major drawbacks: 

1) For each parameter considered the order of the system of equations to 
be solved increases by n. Kahne showed that this order could be somewhat 
reduced by using the fact that the feedback matrix obtained in synnietric. 
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2) The more Important problem Is that in a closed loop Implementation, as 
is usually the case, o(t) no longer represents the trajectory sensitivity. 
This was pointed out by Kreindler [3]. 

We will now formulate the problem for a closed loop implementation. In 
this case u and B will depend on the parameter a. Then eq. (1) gives 

a * X + Ac + B^u + B o(o) * o (B) 

Using the linear feedback law 


we get 


u * K^(t)x + K 2 {t)a 


3u 

3a 


K^o + 



Substituting in eq. (6), 


(7) 

( 8 ) 


6 * (A^ + *^1^^ + (A + Kg + BK^)c + BKg U , ff(o) * 0 (9) 

The second order derivative is neglected because the only way it can 

dOl 

be obtained is by differentiating eq. (9) with respect to a, this introduces 
term and so on. The solution of eq. (9) with neglected will provide 
an approximate trajectory sensitivity vector p(t). 

If p(t) is used in the feedback law, then u » x + KgP (10) 

and substituting in eq. (1) we get 

X = (A + BK,)x + BKgP, x(o) = x^ (11) 

Differentiating (11) with respect to a (neglecting and letting * P) 

gives 

p = (A^ + B^K,)x + (A + BK^ + B^Kg^P* P<0) “ 0 
Exact differentiation of (11) yields 
a = (A^ + B^K,)x + B^KgP + (A + BK^)o + BKg |^ , p(o) * 0 (13) 
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(14) 


where can be obtained from 


li. 

da 


* “a “l’” MA + W, + Kj) ||. || (o) - 0 


The difference between p and o can be reduced by taking u * 1^1 x h- 
where and %£ chosen so as to reduce the difference between p and a. 

Thent 


X ■ (A + ) X + P 

P * (A^ + B^ k^)x + (A + Bk^ + k 2 )p. p(o) • o 

*"» (A„ * B, J(,)a . (A . Bk, * B,^) |t. |fl2l * 0 


(15) 

(16) 
(17) 


A number of authors have used this formulation to solve the problem of 
closed loop trajectory sensitivity. Some of the Important works will now 
be d1 scribed. 

Kreindler [3] used the approximate trajectory sensitivity vector p(t) 
in his formulation. Using equations (1) and (16) he constructed an augmented 
system of order 2n. 


Take 2 *! 
then 


[:] 

t = 5Z + fu, Z(o) = 2 


where 


r A 0 

7 


[(A, + B^ s:,) (A . bJ:, . B^ ^ 

)J * 

[:] 


and 1 


(18) 


also u « [k^ kgjZ 

The augmented problem solved by Kreindler was to find k^ and k 2 so as to 


minimize 


J * r (2^ lyz + u'*’Ru)dt 


(19) 
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subject to (18) 

Where 

Q end S ere positive semldeflnlte nxn metrices end R Is e positive 
definite mxm metrix. 

This Is noM e stenderd llneer reguletor problem end the solution Is given by: 

K « [K^ Kg] “ (20) 

where Is the solution of the metrix Riccetl eq. 

-sj . + lyr - fi^ + ?, i^(t) • o (2D 

As en Initlel estimete for end Kreindler tekes them to be end Kg, 
but then eq. (21) Is no longer one of the Riccetl type beceuse IT depends on K^. 
Kcvertheless he chooses K<| end Kg to setisfy eqs. (20) end (21). 

This eppr^ech works well es long es p does not differ much from the true 
trejectory sensitivity vector a. 

Reo and Soudeck [4] nndifled the above procedure so that one can get a 
better approximation for a. They also use p In the cost functional but use a 
generalized version of Eq.(16). 

p ■ Cx + Dp, p(o) « 0 (22) 

Metrices C end 0 contain as many parameters as one may consider necessary 
to get a close approximation for a by the procedure described below. From eq. 

(22) we get 

” If- IS • 0 

In this case the augmented system of equations are 

Z « AZ + fu, Z(o) » Zq (24) 
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where 


>■[: -[:1 


The procedure Is to choose a set of values for the parameters In C and D 
and with the new expression for X to find a K by solving eqns. (20) and (21). 
This will minimize 0 given by (19). With this K* the parameters In C and D are 


varied so as to minimize 




For this* eqns (11). (13). (22) and (23). have to be Integrated simultaneously 
for a « <>nom1nar values of C and D are used to again solve the Riccatl 

equation to get a new K. This procedure Is repeated until J 2 reduces to an 
acceptable level. 

This procedure may give better results than Kreindler's procedure, but 
Its computational complexity Is far greater than that of the latter. 

Fleming and Newmann [5] use the exact trajectory sensitivity vector 0 
In the performance Index but use the approximate trajectory sensitivity vector p 
In the feedback law. They define the augmented state vector Z as 


P 

Z • 0 

_3- 


Using eqns. (11), (12), (13) and (14) the augmented eqns can be written 


as (of order 4n) 


Z » Z(o) = Z, 


where 


A + BK^ BK 2 0 0 

A * Aqi B(j A + BKi Ba 0 

Aa + Bo Bo K2 A + BK^ BK2 

0 0 * \ h 


1 


and 

0 


0 

0 

0 


Tha performance Index Is 


J (x^Qx + o^sc + u^Ru)dt 


This can be written as 


J •/ (!I + RKE)Z dt 


(26) 


where 


fQ 0 0 0' 
0 0 0 0 
0 0 3 0 
0 0 0 0 


• ^ • [^2n °2n. ^ 


The problem then becomes to choose K so as to minimize J subject to (25). 
Clearly this Is not a linear regulator problem In Its standard form and the 
authors describe a few computational methods of solving It. 

In the three approaches given till now« one either neglects the second 
order derivative term or uses successive minimization techniques. In 

od 

an Interesting correspondence Byrne and Burke [6] use a slightly different 
approach. They use the optimal control for the regulator problem without 
sensitivity constraints as an Initial approximation to the optimal control with 
sensitivity constraints. Let the optimal control for the problem without 
trajectory constraints be 

u * Kx (27) 

where K Is given by 

K«-R"^B^P, (28) 
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( 29 ) 


ind B U the solution of the Riccetl eq. 

P^+ + A^P^ + Q - PjBR*’ B^P^ ■ 0, P(T) ■ 0 

Differentiating u with respect to parameter a gives 

|a-Ka*K„x (30) 

Where K ■ aK/da 
01 

and differentiating (28) with respect to a gives 

K ■ P. (31) 

aP, 

"•“re Z ■ aj 

Differentiating (29) with respect to a 

2+ Z(A-BR‘^B^R|) + (A-BR"^B^Rj)^2 + (P(\^ + A^ R,) 

-P^(B^ R'^B^ + BR“^b 2^) P^ * 0. 2(o) « 0 (32) 

This equation Is solved to g1ve]C> Substituting eq. (30) In eq. (6), we get 

j » (A + BK) a + (A^ + BK^) x + B^u, a(o) * 0 (33) 

Then using the augmented vector Z ■ , the augmented system becomes 

2 » + III, Z(o) * Z^ (34) 

where 


’A 0 


'B * 




. F- 


• ^0 • 


A «► BK A + BK 
a o 




0 

m m 


The problem Is now to minimize the performance Index given by (19) subject 
to (34). This Is a standard regulator problem and the solution Is 

u “ (35) 
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( 36 ) 


Minrt R Is the solution of the equation 

F, ♦ F,X + f F,+ 5 - F^f r' f Fj- 0. F(t) • 0 

A number of comparisons have been made between these methods and the 
results of these depend very much on the problem at hand. But If one was to 
use a comparison Index of some sort composed of reduction In sensitivity per 
unit Increase In cost and the numerical effort required, the.t the routine of 
Kreindler and Byrne and Burke would probably come out ahead. 

Combined trajectory and performance Index sensitivity reduction : The 
parameter sensitivity of the performance Index Is Important because the optimal 
control Is obtained by minimizing the given performance Index. Yahagi [7] 
worked on this problem and gave necessary conditions for an optimal output 
feedback control with reduced performance Index sensitivity. 

Yedavalll and Skelton [1] established relationship between the trajectory 
and performance Index sensitivities. They exploit this relationship to present 
a unified way of reducing trajectory (output or state) sensitivity and performance 
Index sensitivity. Their method also considers control sensitivity, something 
which Is Ignored by most authors. Although Yedavalll and Skelton developed their 
method for the general case of r parameters, we will describe It for the case 
of one parameter only. This somewhat reduces the notatlonal complexity. 

Using the notation developed earlier, the linear time Invariant system Is 


X * Ax ♦ Bu, x(o) * Xq 

(37) 

y ■ Cx 

(36) 


where 

y Is a k dimensional output vector and C Is the output matrix. 
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The performance Index of the unperturbed problem Is 


■f. 


J * |(y^Qy + u^Ru)dt 


(39) 


where 

Q and R are synmetric positive definite matrices. The standard solution of 
this Is 

u * Kx (40) 

where K s - R'^B^R, (41) 

and ^ Is the unique positive definite solution of the matrix Riccatl equation 

Py\ + A^P, - P^R‘^B^P^ + C^QC = 0 (42) 

Let a be the uncertain parameter and A, B, G, x and y are continuous functions 

of a. The sensitivity vectors are 

State trajectory sensitivity vector a * 3x/3a 
Output trajectory sensitivity vector y^ = 

Control trajectory sensitivity vector 
and performance Index sensitivity 

From equations (37) and (38) we get 

a * A X + Aa + B u + Bu^» a(o) = o 
a a a* ' ' 

= C^x + to 


(43) 

(44) 


Mhere the subscript a denotes partlel derivative with respect to a 


Using Z = P and W = ^ the augmented system is 


L<^J L/qiJ 


Z = Az + Su + & u^, Z(o) = Zq 

(45) 

W * Cz 

(46) 
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where 


A 


*A 0' 


'B ■ 


’O' 


. 6 » 


, « 


A 

L (X J 




B 

m* m 


’C 0’ 


m M 

x 



0 


and Z, • 


c c 

•■a J 


0 

l» m 


The authors define a new performance Index which In addition to the output 

sensitivity term also contains control sensitivity term. 

00 

Jg “J {y^Qy + yjQiyci + dt (47) 

where and R-j are positive definite matrices. 

The norms yj Q^ya* '^1^1 functions of time. They showed 

that Is an upper bound to these three norms, 
where 

Is a nonzero scalar such that the matrix 


S 


c 




Is positive definite 


with S * block diagonal [Q, R] 

5 = block diagonal 

2 2 

They also showed that If and are zero, then, 

0 

This relates the trajectory sensitivity terms y' Q, y and uj R,u In (47) with 

w 1 01 Gt I d 

the performance Index sensitivity J^. 

All this Implies that minimization of achieves a trade-off between 
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unperturbed optimal cost J and the sensitivities (output, control and performance 
Index). Although Is not Included In the performance Index J^, performance 
Index sensitivity reduction Is obtained as a "byproduct". 

Using equation (38) and (44) 

Qy + yj Piy^ * x**'c^qcx + (x'*’c^ + aV)Q (c^x + ca) 


then becomes 


■/ 




(C% + cJQiC.) cJq,C 


L c\c 

k- la 


c^q,c 




R 0 


0 R. 


11 L^oJ) 


Idt 


(49) 


The above expression requires u^. Since the desired control law 
u = k-jx + k 2 d 


(50) 


Is not available, the authors use the control for the unperturbed problem 
(given by eqs. (40), (41) and (42)) to obtain an approximate u^. 

Then the approximate control sensitivity u^ is given by 

% • kP (51 ) 

where p is the approximate trajectory sensitivity. 

Substituting equation (51) in eq. (43), 

p » A^x + (A + BK)p + B^u (52) 


Defining Z. * 

a 


, we get 


u. = z. 


a. 


(53) 

(54) 
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where 



o 

t 


’B ■ 


X 

0 

*1 ■ 


f » 


• h • 



A, (A + BK). 


•“a- 

0 

.0 . 


and ^ * [0 K] 


The performance Index J. becomes (using for u ) 

s Ola 01 


f‘ 


Js • I h * « 


(55) 


where U, 


(c^Qc + c^Oi g cj(^c 

Cjj (c^Q, c + Jr, k). 


The problem Is to minimize subject to equation (53). This Is a standard 
regulator problem and the solution Is 

u = RZ^*RjX + K2P (56) 

where R * R| (57) 

fi| is the positive definite solution of the Riccatl equation 

P^l + a||^ - = 0 (58) 

So the steps in the Yedavalli and Skelton procedure are 

(a) Compute K given by eqns. (41) and (42) 

(b) Form the matrices A^, and 

(c) The desired control is given by equations (56), (57) and (58). 

As pointed earlier, results of comparisons between the various methods 
depends on the problem used. Keeping this In mind we will now use a first 
order example for the comparison of the methods described above. This example 
has been used by a number of authors for this purpose. 
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The first order plant is 
, 2 

X = O X + u, x(o) = 1 

where a is the parameter with a nominal value of 1. The objective is to minimize 
the trajectory sensitivity. 

The performance index is given by 

J = J(x^ + u^) dt 
0 

To aid in the comparison of the various methods, the sensitivity integral, 

S, is also evaluated along with J for each of the methods. S is given by 



The results of this are given in table 1. Clearly for all the algorithms 
a comparison with the simple regulator reveals that the decrease in sensitivity 
measured by S is obtained at the expense of a higher cost given by J. The best 
method would be the one which for the least increase in J will give the maximum 
decrease in S. The method of Byrne and Burke is the best for this example. 

The methods of Rao and Saudack and Fleming and Newman require much more 
computation than the other methods. 
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simple Regulator Rao and Fleming and Byrne and Yedavalll and 
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Eigenvalue Sensitivity 

Besides the trajectory and performance Index (cost) sensitivities there 
Is another Important means of estimating the sensitivity of the system to 
variations In Its parameters. This Is through the use of eigenvalue/eigenvector 
sensitivity. It should be mentioned here that eigenvalue/eigenvector sensitivity 
provides a less direct measure of the system sensitivity than trajectory 
or cost sensitivity. Designers working with classical design techniques tend 
to use eigenvalue sensitivity. 

Crossby and Porter [8] developed expressions for eigenvalue and eigenvector 
sensitivities for linear time Invariant systems. Reddy [9] also worked on 
the problem of determining the effects of variation In system parameters on 
eigenvalues. Although they provided explicit expressions for the sensitivity, 
these authors do not suggest a method of reducing the eigenvalue/eigenvector 
sensitivity. In an Interesting paper Gourishankar and Ramar [10] combine the 
problems of reducing eigenvalue sensitivity and closed loop eigenvalues (pole) 
placements for a linear time Invariant multivariable system. We will now discuss 
their approach In some detail. 

It has been known for a long time that using complete state feedback 
the closed loop poles can be assigned to any desired location. It Is also 
known that for multi Input systems a number of feedback matrices give the desired 
closed loop pole locations. These controllers will In general give different 
time responses. Designers usually use this freedom to obtain other desired 
characteristics in the system response. Gourishankar and Ramar use this 
freedom to minimize the sensitivity of the closed loop eigenvalues to the variation 
In the system parameters without effecting their desired location. In the 
following development, the closed-loop poles are taken to be distinct. 
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Using the notation of the previous chapter, the linear time Invariant 
system Is given by 

X * Ax + Bu (1) 

and the feedback law Is 

u » kx (2) 

where K Is mxn time Invariant feedback matrix. 

The closed loop system Is then given by 

X « (A + BK)x (3) 

Keeping In mind that n elements of the mxn feedback matrix are sufficient 
to place the n closed loop poles to their desired location, we subdivide 
the feedback matrix Into two parts as follows 

R » [kj kj k^.j k^ J kj]"^ 

/s 

and k = k^. 

where kj Is the jth row of the K matrix, j » 1, 2, — m 

A 

The vector k Is used to obtain the desired pole locations while the 
elements of the matrix R are chosen to reduce the eigenvalue sensitivity. 
Define 

5 • [u, Uj 

A 

and u * u^. 

Eq. (2) becomes 

u * Rx (4) 

and 

u » kx (5) 
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We can now write eq. (1) as 

X ■ (A + B R)x + b u (6) 

where 

6 - [b^ bg b^., b^^^ b^] and b » b^ 

bj's being the columns of the matrix B. 

A. 

Eq. (6) represents a single Input system (u « u^} and the feedback vector 
k In eq. (5) can be designed to place the closed loop poles In their specified 

A 

location. This can always be done as long as the pair (A ■»' § K, b) Is completely 
controllable. Since A -t* 6 R depends on R, care has to be taken while selecting 
R so as to maintain the above mentioned pair completely controllable. The 
authors have mentioned that this is not a serious limitation as almost any R 
satisfies this requirement. 

Substituting eq. (5) in (6} yields the closed loop system 
X = (A + § R + b k)x (7) 

A 

It is to be noted that the feedback vector k depends on A + 6 K and hence 

on K. 

Morgan [11] gave the expression for the sensitivity of the closed loop 
eigenvalues to the variations in the elements of the system matrix A. 

■ TlgTsTTisF^ 

where 

Sjl sensitivity of the eigenvalue to a small variation in the element 

of the system matrix A, 

g(s) is the characteristic polynamial of the closed loop system, 

R(S^)= adjoint (S^ I - %) 

A A 

and A*A + BR + bk 
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The authors choose the elenients of the matrix R so as to minimize the 


performance index 


J 


n n n 

I I I 




(9) 


If only the a^^ element is changing, the performance index becomes 

J • I (10) 

1*1 

The authors mention that transforming to phase- variable form facilitates in 
the minimization of J. 

The procedure can be summarized as follows 

(a) Choose R so as to minimize J 

A 

(b) With this value of R find k which assigns the closed loop poles to the 
desired location. 

A 

Any row of matrix R can be used as k. One could take the rows of R as 

A 

k one at a time and calculate J^^^imum choose that row of R 

A 

as k for which the lowest was obtained. This is usually not done as 

the amount of computation increases enormously 

Tne design procedure is illustrated by tie following linear time invariant 
second-order system with two inputs and two outputs, 

X * Ax + Bu (11 ) 

where 



1 

p* 

ro 

1 


'1 

r 

A » 


. B = 




.*21 * 22 - 


0 

1 , 


The nominal values of a-|^, ai 2 > a 2 | and a 2 £ 0, 1, 0 and 0 respectively. 
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Only 82^ varies from Its nominal value. The desired location for the closed 
loop poles is » -2 and $2 * - 3 . The feedback matrix Is 


K 


A 


hi 

1 

CM 

j 2 “ 

X 

•k* 

•hi 

1 

CM 

CM 


.R 


then 


B • [b B] 


and % 


This gives 


■ 0 r 

+ 

■ 1 * 

[k2i k22l + 

» « 
1 

.*21 


1 


0 

■ m 

■(K„ * 

•^21 

) 

(1 + k-|2 + k2£) 

.(»21 * 

k2i 

) 

CM 

CM 



'0 

O' 



^*21 

1 

m 

0, 





r(s 

- k^^ - k2^) 

-i 


Ck^l k^2^ 


and R(s) * Adjoint 


* ^*21 ^ 

* ^ 22 ^ ^ ^12 ^ 22 ^ 

^*21 “ hi “ hi^ 


T ,.^2 - -^22^ 

(s - k22) 
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Now trace R(s) 


3% 


trt s 


(1 ^ ^12 * ^ 22 ^ 
(s - . kg^) 


0 ’ 

0 


• (1 ♦ k ^2 ■** -^ 22 ^ 


9 (s) • (s ■ k-|i - ^21^ " ^22^ • (1 + k -|2 ’*’ ^22^ ^*21 ^ ^21^ 


.'. ■ 2s - (kll ♦ 1^1 + •‘22) 


The eq. (8) now gives 




-2 


"(1 * k ^2 ■*■ ^ 22 ^ 

(k^^ + k 2 ^ + k 22 + 4 ) 



as 



-3 


«(1 •«■ k ^2 ^22^ 

(k|^ + k 2 | + k 22 + 6 ) 


The performance Index becomes 


J • (S 2] )2 + 

* (1 + k *2 ^ ^ 2 ^ I ■ 2 ^ " ' " 2 

l(k^^ + k2i + k22 '*’ 4) (k^l + k2-j + k2g + 6) 


We have to choose R so as to minimize J. J takes Its minimum value 0 If 


k22 ■ "(1 *► k^2^ 

k^l + k2i + k22 ♦ 4 f 0 
and k^^ k 2 ^ + k 22 + 6 f 0 
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Take 1^] * 0 
For * 0, we get 

g($) »“ 1^21 *^22^ *^n *^22 “ *^21 * Hz 4l 

Substituting the values of k^^ and k 2 £ In the above expression* 
g(s) • s^ -s(k,^ - k ^2 ‘ * *^12^ 

k^l* k ^2 ^ ^ chosen so as to get the desired characteristic polynomial 

(s + 2) (s + 3) * s^ + 5s + 6 
Equating the two polynomials, we get 

- 1^12 * 1 ■ ‘5 

and k^^ (1 + k^ 2 ) * 

Solving for k^ 2 * 

(k,2 - 4) (1 + k,2> * 
or k ^2 *3 k ^2 ■*“ 2 ■ 0 

or (k ^2 ‘ 2) (k ^2 * 1) * 0 

k„ . -2. -3 
k22 * *3, “2 
‘‘21 ■« 

This gives the desired feedback matrix as 



1 

CM 

CM 

1 

1 


■-3 r 

K * 


or 



0 -3 

m ^ 


1 

0 

1 

ro 
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ROBUST CONTROLLER DESIGN FOR LARGE PARAMETER VARIATIONS 
Leonard R. Anderson* 


ASSUMPTIONS 

1) The physical process depends on a few key parameters ay ...» oj^. 
(For wing flutter suppression these would be dynamic pressure, fuel 
load In wings, etc). N would not typically be larger than 2 or 3. 

2) These N parameters vary c er limits that are either specified , or to 
be determined in the design process. 


ttj - i * 1 , . . . , N 

‘min ‘ ‘max 

The range of each parameter Is divided into a grid of a small number 
of equally spaced points, e.g. , 3-7 points. 


[oj ; a. ] = [a. ; a. *, 

'nin max 1 2 



] 


with a. = aj ; o. = a. 

‘min ‘l ‘max ‘Mi 

The set of of possible parameter values, therefore has N^ * M^ • Mj, 
elements. 

3) For each element ^ e the physical process can be modeled by a system 
of linear constant coefficient O.O.E.'s. For example, consider a wing 
flutter model. 


♦Assist. Prof., Aerospace and Ocean Eng. Dept., VPI & SU. 
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X • A(t)x + B(a)u G(t) w. teR". u R™ 
y * Hx + V yeR*^ 



Z = Ac? + B^y 
iT = K^? + K2? 


Typical Parameters that would vary might be 

Actuators: Gain, Phase Error 
Vehicle : Dynamic Pressure, Geometry 
Sensors : Gain, Offset Error 


4) The structure of the control law is specified, but values of constants 
in the control law are to be determined in the design process. 

For example, suppose the sensor signals consist of accelerations, y^, 
measured at various points on the vehicle/wing, and the control law is: 



u = t + K 2 y 

Let the elements of matrices K 2 and C-|-j, € 2 ^* ^22* “ 

represented by the design vector iT. 

5) If we choose a parameter vector a and a design vector iT, we have a 
well defined closed-loop linear system 
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6) For one design point k there are closed- loop system matrices 

^ -fc ''' j» 

A(k$oi<)j •••» A(k$oini ) 

o 

corresponding to the possible discrete values of the parameters 
* ^2 * • * • • 

a 

We can compute the eigenvalues of each of these closed-loop system 
matrices as 

A (A(M^)) 

and collect them all Into one set that depends only on the design 
point k as 

A (k) = |A{A(k,?,)) hI(X,Z. ))( 

'a 

7) Also, for specified Initial state and noise covariance matrices, we 
can compute the RMS values of actuator signals, structural degrees 

of freedom, or other linearly related parameters for any given closed 
loop system A(k,a). 

If the number were not unreasonably large, we could find the maxi mum 

RMS values of these variables for the closed-loop systems 

^ 

A(kjOin)* A(k»cxjj ) 

8) Then, applying either random pattern search strategy (Bekey *81) or 
nonlinear programming without gradient evaluation, we can search over 
the space S]^of possible design points to achieve any of the following 
design objectives: 

a) minimize RMS values of actuator signals while requiring the set 
A(iT) to remain In some favorable region of the complex plane. 

b) maximize the range while constraining both RMS 

values of actuator signals, and location of A(k) In the complex plane. 
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j 

I c) for specified limits on actuator signals and parameter ranges, move 

a(|T) as far to the left in the complex plane as possible. 

IMPLEMENTATION PROBLEMS; 

1) One must be able to effectively generate the linear models A(^, 

B®, G(a) corresponding to the specified parameter ranges. This is 
probably not possible on-line, but will require these arrays of arrays 
be stored on disk before the controller design study is initiated. 

2) It may be difficult to avoid local minima in the search for an optimal 
Ic, and alternate starting points should be considered. 

3) One must religiously avoid the curse of dimensionality and reduce the 
number of parameter points N^ to an absolute minimum (i.e., 9 or 15). 

4) Care must be taken to construct a well-posed problem in which none of 
the design variables t may go to infinity. 

REF; G. A. Bekey, "Modeling and Identification of Nonlinear Systems 

( 

with Hysteresis", 3rd Int. Conf. on Math. Modeling, U.C.L.A., July 
29-31 , 1981 . 
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